library(tidyverse)
library(ISLR)
library(kableExtra)
library(GGally)
library(knitr)
library(lmtest)
library(corrplot)
library(olsrr)
library(caret)
library(plotly)
house <- read.csv('C:/Users/Kuba/Downloads/house.csv', header = TRUE)
head(house)
Zbiór house przedstawia ceny domów w Stanach
Zjednoczonych. Niektóre istotne kolumny w zbiorze danych:
AvgAreaIncome - średni dochód powierzchniowy
mieszkańców domu w [$],HouseAge - wiek domu,NumberOfRooms - liczba pokoi (średnia),NumberOfBedooms - liczba sypialni (średnia),AreaPopulation - liczba ludności w dzielncy/regionie, w
której znajduje się dany dom,Price - cena domu w [$].kable(summary(house), escape = F) %>%
kable_styling(bootstrap_options = c("striped", "hover", "condensed"))
| Avg..Area.Income | House.Age | Number.of.Rooms | Number.of.Bedrooms | Area.Population | Price | Address | |
|---|---|---|---|---|---|---|---|
| Min. : 17797 | Min. :2.644 | Min. : 3.236 | Min. :2.000 | Min. : 172.6 | Min. : 15939 | Length:4548 | |
| 1st Qu.: 61485 | 1st Qu.:5.332 | 1st Qu.: 6.300 | 1st Qu.:3.140 | 1st Qu.:29423.2 | 1st Qu.: 997775 | Class :character | |
| Median : 68817 | Median :5.961 | Median : 7.002 | Median :4.050 | Median :36215.6 | Median :1234571 | Mode :character | |
| Mean : 68612 | Mean :5.979 | Mean : 6.988 | Mean :3.982 | Mean :36187.5 | Mean :1233916 | NA | |
| 3rd Qu.: 75821 | 3rd Qu.:6.658 | 3rd Qu.: 7.666 | 3rd Qu.:4.490 | 3rd Qu.:42880.6 | 3rd Qu.:1470616 | NA | |
| Max. :107702 | Max. :9.519 | Max. :10.760 | Max. :6.500 | Max. :69592.0 | Max. :2469066 | NA |
Przekształcamy “tibbl’a” ze wszystkimi zmiennymi na mniejszą
wersję, zawierającą wyłącznie zmienne ciągłe (usuwamy kolumnę
Address).
data(house)
## Warning in data(house): data set 'house' not found
ggpairs(house[, c(1:6)])
Powyższy wykres przedstawia zależności pomiędzy wszystkimi ciągłymi
zmiennymi. Można zauważyć, że istnieje kilka znacznych korelacji. Jedna,
która wyróżnia się ponad innymi to ta pomiędzy zmienną
Price i AvgAreaIncome ze współczynnikiem
korelacji Pearsona równym \(0,641\).
Korelacje te są dość istotne, ponieważ wartość współczynnika jest bliska
\(1\) dla współczynnika dodatniego.
Dlatego ceny domów w Stanach Zjednoczonych będziemy przewidywać na
podstawie zmiennej Avg..Area.Income. Innymi słowy zmienna
Price to zmienna objaśniana, a
Avg..Area.Income objaśniająca.
Sprawdzimy teraz rozkład tych dwóch zmiennych za pomocą histogramów.
ggplot(house, aes(x=Price))+
labs(title = 'Histogram cen domów',
x = 'Cena',
y = 'Częstotliwość') +
geom_histogram(aes(y=..density..),
color="khaki1",
fill="mediumslateblue",
binwidth = 50000) +
geom_density(alpha=.45, fill="springgreen4") +
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## i Please use `after_stat(density)` instead.
ggplot(house, aes(x=Avg..Area.Income))+
labs(title = 'Histogram średniego dochodu powierzchniowego w domu',
x = 'Średni dochód powierzchniowy w domu',
y = 'Częstotliwość') +
geom_histogram(aes(y=..density..),
color="khaki1",
fill="mediumslateblue",
binwidth = 2100) +
geom_density(alpha=.45, fill="springgreen4") +
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
Wizualizując rozkłady tych dwóch zmiennych możemy domyśleć się, że posiadają rozkład normalny. Dla potwierdzenia naszej hipotezy użyjemy testu Shapiro-Wilka, który zakłada, że zmienna posiada rozkład normalny.
shapiro.test(house$Price)
##
## Shapiro-Wilk normality test
##
## data: house$Price
## W = 0.9998, p-value = 0.9688
shapiro.test(house$Avg..Area.Income)
##
## Shapiro-Wilk normality test
##
## data: house$Avg..Area.Income
## W = 0.99958, p-value = 0.4456
Wnioskując, w obu przypadkach otrzymaliśmy \(p\)-value \(> 0,05\), zatem nie mamy podstaw do odrzucenia \(H_0\).
scatter_plot <- plot_ly(data = house, x = ~Avg..Area.Income, y = ~Price, color = ~Price) %>%
layout(title = 'Wykres punktowy zależności ceny domu od średniego dochodu powierzchniowego', plot_bgcolor = "#e5ecf6")
scatter_plot
## No trace type specified:
## Based on info supplied, a 'scatter' trace seems appropriate.
## Read more about this trace type -> https://plotly.com/r/reference/#scatter
## No scatter mode specifed:
## Setting the mode to markers
## Read more about this attribute -> https://plotly.com/r/reference/#scatter-mode
Wykres nie wyklucza zależności liniowej pomiędzy zmiennymi
Price i Avg..Area.Income, więc przejdziemy do
stworzenia modelu regresji liniowej tych zmiennych.
lmfit1 <- lm(Price~Avg..Area.Income, data = house)
summary(lmfit1)
##
## Call:
## lm(formula = Price ~ Avg..Area.Income, data = house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -965631 -187471 -3066 184157 982367
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.262e+05 2.622e+04 -8.627 <2e-16 ***
## Avg..Area.Income 2.128e+01 3.775e-01 56.365 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 272100 on 4546 degrees of freedom
## Multiple R-squared: 0.4114, Adjusted R-squared: 0.4112
## F-statistic: 3177 on 1 and 4546 DF, p-value: < 2.2e-16
Z naszego modelu, możemy wywnioskować:
ggplot(lmfit1, aes(x=resid(lmfit1)))+
labs(title = 'Histogram reszt modelu',
x = 'Reszty',
y = 'Częstotliwość') +
geom_histogram(aes(y=..density..),
color="khaki1",
fill="mediumslateblue",
binwidth = 50000) +
geom_density(alpha=.45, fill="springgreen4") +
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
ggplot(lmfit1, aes(sample=lmfit1$residuals)) +
geom_qq() +
geom_qq_line(color = 'red') +
labs(title='Wykres kwartyl-kwartyl reszt', x='Kwartyle teoretyczne', y='Kwartyle próbkowe')+
theme(plot.title = element_text(hjust = 0.5))
shapiro.test(lmfit1$residuals)
##
## Shapiro-Wilk normality test
##
## data: lmfit1$residuals
## W = 0.99971, p-value = 0.7937
Wszystkie metody użyte przez nas świadczą o tym, że reszty naszego modelu posiadają rozkład normalny.
Do sprawdzenia zerowej średniej reszt użyjemy testu t-studenta.
t.test(lmfit1$residuals)
##
## One Sample t-test
##
## data: lmfit1$residuals
## t = -9.1167e-15, df = 4547, p-value = 1
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
## -7908.126 7908.126
## sample estimates:
## mean of x
## -3.677457e-11
W tym wypadku test t-student wykazał, że średnia jest równa zero.
ggplot(lmfit1, aes(.fitted, .resid)) +
geom_point() +
stat_smooth(method='loess', formula=y~x, se=F) +
geom_hline(yintercept=0, linetype='dashed', color='red') +
labs(title='Wykres zależności reszt od dopasowanych wartości',
x='Dopasowane wartości',
y='Reszty')+
theme(plot.title = element_text(hjust = 0.5))
Patrząc na wykres zależności reszt od dopasowanych wartości możemy
zauważyć większe odchylenia dla małych jak i dużych wartości zmiennej
objaśnianej (Price).
W celu sprawdzenia niezależności reszt przeprowadzimy test Durbin-Watson’a.
lmtest::dwtest(lmfit1)
##
## Durbin-Watson test
##
## data: lmfit1
## DW = 1.9676, p-value = 0.1373
## alternative hypothesis: true autocorrelation is greater than 0
W naszym przypadku \(p\)-value \(> 0,05\), więc nie mamy podstaw, aby odrzucić hipotezę o niezależności w resztach.
Aby sprawdzić homoskedastyczność posłuży nam wykres zależności pierwiastka standaryzowanych reszt od dopasowanych wartości.
lmfit1 %>%
ggplot(aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point() + stat_smooth(method='loess', formula=y~x, se=F) +
labs(title='Zależność pierwiastka standaryzowanych reszt od dopasowanych wartości', x='Dopasowane wartości', y='Pierwiastek standaryzowanych reszt')+
theme(plot.title = element_text(hjust = 0.5))
Rozrzut reszt na wykresie jest mniej więcej równy dla wszystkich dopasowanych wartości. Widzimy odchylenie na początku niebieskiej prostej. Jest ono spowodowane wartością odstającą. Deszty są równomiernie rozrzucone wokół niebieskiej linii. Dla pewności przeprowadzimy test Breusch-Pagan’a, który za hipotezę zerową zakłada homoskedastyczność reszt.
lmtest::bptest(lmfit1)
##
## studentized Breusch-Pagan test
##
## data: lmfit1
## BP = 1.893, df = 1, p-value = 0.1689
\(P\)-value \(> 0,05\), więc zakładamy homoskedastyczność reszt.
Słowem podsumowania stwierdzamy, że powyższy model regresji liniowej
jest zgodny z naszymi założeniami. Korzystająć ze współczynnika
korelacji Pearsona dopasowaliśmy zmienną objaśniającą
(Avg..Area.Income). Przewidywana zmienna Price
jest najbardziej zależna od zmiennej Avg..Area.Income.
Analiza reszt pokazała, że model jest zgodny założeniami i ma
postać:
\(Price = -226200 + 21,28 \cdot Avg..Area.Income\)
Spójrzmy raz jeszcze na nasz zbiór danych.
head(house)
Ze względnu na to, że w naszym modelu nie może byc zmiennych
jakościowych, wybieramy tylko zmienne ilościowe przy pomocy funkcji
select().
new_house <- house %>%
select(Avg..Area.Income, House.Age, Number.of.Rooms, Number.of.Bedrooms, Area.Population, Price)
head(new_house)
lmfit2 <- lm(Price~., data = new_house)
summary(lmfit2)
##
## Call:
## lm(formula = Price ~ ., data = new_house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -338141 -70192 -101 69104 347218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.639e+06 1.791e+04 -147.400 <2e-16 ***
## Avg..Area.Income 2.161e+01 1.404e-01 153.929 <2e-16 ***
## House.Age 1.653e+05 1.514e+03 109.228 <2e-16 ***
## Number.of.Rooms 1.211e+05 1.682e+03 72.003 <2e-16 ***
## Number.of.Bedrooms 1.458e+03 1.377e+03 1.059 0.289
## Area.Population 1.519e+01 1.513e-01 100.357 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 101100 on 4542 degrees of freedom
## Multiple R-squared: 0.9188, Adjusted R-squared: 0.9187
## F-statistic: 1.028e+04 on 5 and 4542 DF, p-value: < 2.2e-16
Dzięki funkcji summary() odczytujemy, że \(R^2\) naszego modelu wynosi \(0,9187\). Tak zbudowany model wyjaśnia
\(91,87\%\) zmienności bieżącej ceny,
ale nie wszystkie zmienne są w tym modelu istotne.
Zinterpretujmy nasze wyniki:
Wykorzystamy współczynnik informacyjny Akaike'a
(AIC). Jest to estymator błędu przewidywania. Biorąc pod
uwagę zbiór modeli dla danych, AIC szacuje jakość każdego
modelu, w stosunku do każdego z pozostałych modeli. W ten sposób
AIC dostarcza środków do wyboru modelu. Dzięki niemu
uzyskamy najlepszą możliwą dla nas jakość modelu poprzez eliminację
następnych zmiennych. Skorzystamy z powyższego algorytmu przy pomocy
funkcji step(). Algorytm będzie działał na zasadzie
procedury eliminacji wstecznej.
step(lmfit2, direction = "backward")
## Start: AIC=104827.4
## Price ~ Avg..Area.Income + House.Age + Number.of.Rooms + Number.of.Bedrooms +
## Area.Population
##
## Df Sum of Sq RSS AIC
## - Number.of.Bedrooms 1 1.1473e+10 4.6440e+13 104827
## <none> 4.6429e+13 104827
## - Number.of.Rooms 1 5.2996e+13 9.9425e+13 108289
## - Area.Population 1 1.0295e+14 1.4938e+14 110140
## - House.Age 1 1.2196e+14 1.6839e+14 110685
## - Avg..Area.Income 1 2.4220e+14 2.8863e+14 113136
##
## Step: AIC=104826.6
## Price ~ Avg..Area.Income + House.Age + Number.of.Rooms + Area.Population
##
## Df Sum of Sq RSS AIC
## <none> 4.6440e+13 104827
## - Number.of.Rooms 1 6.8538e+13 1.1498e+14 108948
## - Area.Population 1 1.0295e+14 1.4939e+14 110138
## - House.Age 1 1.2204e+14 1.6848e+14 110685
## - Avg..Area.Income 1 2.4255e+14 2.8899e+14 113139
##
## Call:
## lm(formula = Price ~ Avg..Area.Income + House.Age + Number.of.Rooms +
## Area.Population, data = new_house)
##
## Coefficients:
## (Intercept) Avg..Area.Income House.Age Number.of.Rooms
## -2.640e+06 2.161e+01 1.654e+05 1.220e+05
## Area.Population
## 1.519e+01
Jak widzimy kolejno usuwamy zmienne, które w naszym modelu mają
najmniejszą wartość AIC, najpierw był to
Number.of.Bedrooms. Dochodzimy do momentu w algorytmie gdzie
najniższą wartość AIC ma <none> zatem
poniższe zmienne zostają w naszym modelu regresji wielorakiej. Widzimy,
że dla wszystkich zmiennych \(R^2\)
wynosi \(0,9187\). Sprawdźmy teraz ile
wynosi nasz \(R^2\) dla naszego
gotowego modelu, po redukcji tej zmiennej zmiennych:
new_lmfit2 <- lm(Price~ Avg..Area.Income + House.Age + Area.Population + Number.of.Rooms, data= new_house)
summary(new_lmfit2)
##
## Call:
## lm(formula = Price ~ Avg..Area.Income + House.Age + Area.Population +
## Number.of.Rooms, data = new_house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -339350 -70321 -126 69246 347378
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.640e+06 1.790e+04 -147.45 <2e-16 ***
## Avg..Area.Income 2.161e+01 1.403e-01 154.04 <2e-16 ***
## House.Age 1.654e+05 1.513e+03 109.27 <2e-16 ***
## Area.Population 1.519e+01 1.513e-01 100.35 <2e-16 ***
## Number.of.Rooms 1.220e+05 1.490e+03 81.88 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 101100 on 4543 degrees of freedom
## Multiple R-squared: 0.9188, Adjusted R-squared: 0.9187
## F-statistic: 1.284e+04 on 4 and 4543 DF, p-value: < 2.2e-16
Nasz współczynnik \(R^2\) dla obu
modeli jest cały czas taki sam, zatem zmienna
Number.of.Bedrooms nie ma wpływu istotnie na jakość modelu.
Natomiast pozostałe zmienne są istotne dlatego nie odrzuciliśmy ich przy
użyciu algorytmuAIC.
new_lmfit2_no_avg <- lm(Price~ House.Age + Area.Population + Number.of.Rooms, data=new_house)
summary(new_lmfit2_no_avg)
##
## Call:
## lm(formula = Price ~ House.Age + Area.Population + Number.of.Rooms,
## data = new_house)
##
## Residuals:
## Min 1Q Median 3Q Max
## -895558 -167437 2368 171008 878807
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.126e+06 3.732e+04 -30.16 <2e-16 ***
## House.Age 1.630e+05 3.775e+03 43.18 <2e-16 ***
## Area.Population 1.501e+01 3.774e-01 39.77 <2e-16 ***
## Number.of.Rooms 1.205e+05 3.715e+03 32.42 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 252200 on 4544 degrees of freedom
## Multiple R-squared: 0.4945, Adjusted R-squared: 0.4941
## F-statistic: 1481 on 3 and 4544 DF, p-value: < 2.2e-16
Jak widzimy powyżej, nasz \(R^2\) znacznie zmalał, co oznacza iż ten nowy model objaśnia tylko \(49,41\%\) zmienności bieżacej ceny domów. Jest to niekorzystne w wypadku gdy możemy objaśniać aż \(91,87\%\) ceny.
Omówmy temat współliniowości zmiennych objaśniających:
ols_vif_tol(new_lmfit2)
Kolumna Tolerance wskazuje wartość procentową
niewyjaśnionej zmienności danej zmiennej przez pozostałe zmienne
objaśniające. Dla przykładu: współczynnik tolerancji dla
Avg..Area.Income wynosi $0.9997982 $ co oznacza, że \(99,97\%\) zmienności
Avg..Area.Income nie jest wyjaśnione przez pozostałe
zmienne w modelu. Kolumna VIF jest obliczana na podstawie
wartości współczynnika tolerancji i pokazuje o ile wariancja szacowanego
współczynnika regresji jest podwyższona z powodu współliniowości danej
zmiennej objaśniającej z pozostałymi zmiennymi objaśniającymi. Wszystkie
zmienne mają VIF na podobnym poziomie w granicach \(1.0\) zatem nie wskazują współliniowości
(VIF \(> 4\) wskazuje
na współliniowość zmienncyh, także u nas takie zmienne nie
występują).
Sprawdźmy normalność naszych reszt:
ggplot(lmfit2, aes(x=resid(lmfit2)))+
labs(title = 'Histogram reszt modelu',
x = 'Reszty',
y = 'Częstotliwość') +
geom_histogram(aes(y=..density..),
color="khaki1",
fill="mediumslateblue",
binwidth = 18000) +
geom_density(alpha=.2, fill="springgreen4") +
theme_classic()+
theme(plot.title = element_text(hjust = 0.5))
Reszty w naszym modelu wydają się być zbliżone do rozkładu normalnego, natomiast nie możemy tego stwierdzić w 100%. Sprawdźmy wykres QQ-plot do zbadania rozkładu normalnego naszych zmiennych.
ggplot(new_lmfit2, aes(sample=new_lmfit2$residuals)) +
geom_qq() +
geom_qq_line(color = 'red') +
labs(title='Wykres kwartyl-kwartyl reszt', x='Kwartyle teoretyczne', y='Kwartyle próbkowe')+
theme(plot.title = element_text(hjust = 0.5))
Widzimy, że prawie wszystkie punkty leżą na czerwonej lini, jedynie problem jest na początku oraz na końcu.
new_lmfit2 %>%
ggplot(aes(.fitted, .resid)) +
geom_point() +
stat_smooth(method='loess', formula=y~x, se=F) +
geom_hline(yintercept=0, linetype='dashed', color='red') +
labs(title='Wykres zależności reszt od dopasowanych wartości',
x='Dopasowane wartości',
y='Reszty')+
theme(plot.title = element_text(hjust = 0.5))
Na wykresie zależności reszt od dopasowanych wartości możemy zauważyć delikatne odchylenie w początkowej jego fazie. Innymi słowy, współczynniki modelu regresji powinny być godne zaufania i nie musimy wykonywać transformacji na danych.
ols_plot_cooksd_bar(new_lmfit2)
ols_plot_resid_stud_fit(new_lmfit2)
Przypisanie tej funkcji do obiektu zwraca nam tabelę z numerami zidentyfikowanych obserwacji wpływowych. Sposób ten pozwala nam na identyfikację punktów, które negatywnie wpływają na model regresji są one oznaczone kolorem czerwonym. Miara jest kombinacją wartości dźwigni i reszt każdej obserwacji; im wyższa dźwignia i reszty, tym wyższa odległość Cooka.
new_lmfit2 %>%
ggplot(aes(.fitted, sqrt(abs(.stdresid)))) +
geom_point() + stat_smooth(method='loess', formula=y~x, se=F) +
labs(title='Zależność pierwiastka standaryzowanych reszt od dopasowanych wartości', x='Dopasowane wartości', y='Pierwiastek standaryzowanych reszt')
Na powyższym wykresie zależności pierwiastka standaryzowanych reszt od dopasowanych wartości, możemy zauważyć, iż wartości są rozmieszczone mniej więcej w takim samym odstępie od niebieskiej linii. Nie wykluczamy, więc homoskedastyczności. Dla pewności przeprowadźmy test Breusch-Pagan’a:
lmtest::bptest(new_lmfit2)
##
## studentized Breusch-Pagan test
##
## data: new_lmfit2
## BP = 5.8309, df = 4, p-value = 0.2121
Nie mamy wystarczających dowodów, aby stwierdzić, że heteroskedastyczność występuje w naszym modelu regresji, zatem możemy założyć homoskedastyczność.
Stworzony przez nas model spełnia wszystkie założenia. Korzystając z
algorytmu AIC dobraliśmy najlepszy możliwy model regresji
wielorakiej, który dany jest wzorem:
\(Price = -2640000 + 21,61 \cdot Avg..Area.Income + 165400 \cdot House.Age +\\+ 15,19 \cdot Area.Population + 122000 \cdot Number.of.Rooms\)